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Abstract.  The  central  scheme  of  Nessyahu  and  Tadmor  [J.  Comput.  Phys,  87  (1990)]  solves  hyper¬ 
bolic  conservation  laws  on  a  staggered  mesh  and  avoids  solving  Riemann  problems  across  cell  bound¬ 
aries.  To  overcome  the  difficulty  of  excessive  numerical  dissipation  for  small  time  steps,  the  recent  work 
of  Kurganov  and  Tadmor  [J.  Comput.  Phys,  160  (2000)]  employs  a  variable  control  volume,  which  in 
turn  yields  a  semi-discrete  non-staggered  central  scheme.  Another  approach,  which  we  advocate  here, 
is  to  view  the  staggered  meshes  as  a  collection  of  overlapping  cells  and  to  realize  the  computed  solution 
by  its  overlapping  cell  averages.  This  leads  to  a  simple  technique  to  avoid  the  excessive  numerical 
dissipation  for  small  time  steps  [Y.  Liu;  J.  Comput.  Phys,  209  (2005)].  At  the  heart  of  the  proposed 
approach  is  the  evolution  of  two  pieces  of  information  per  cell,  instead  of  one  cell  average  which  charac¬ 
terizes  all  central  and  upwind  Godunov-type  finite  volume  schemes.  Overlapping  cells  lend  themselves 
to  the  development  of  a  central- type  discontinuous  Galerkin  (DG)  method,  following  the  series  of  work 
by  Cockburn  and  Shu  [J.  Comput.  Phys.  141  (1998)]  and  the  references  therein. 

In  this  paper  we  develop  a  central  DG  technique  for  hyperbolic  conservation  laws,  where  we  take 
advantage  of  the  redundant  representation  of  the  solution  on  overlapping  cells.  The  use  of  redundant 
overlapping  cells  opens  new  possibilities,  beyond  those  of  Godunov-type  schemes.  In  particular,  the 
central  DG  is  coupled  with  a  novel  reconstruction  procedure  which  post-processes  the  central  DG 
solution  to  remove  spurious  oscillations  in  the  presence  of  shocks.  This  reconstruction  is  motivated 
by  the  moments  limiter  of  Biswas,  Devine  and  Flaherty  [Appl.  Numer.  Math.  14  (1994)],  but  is 
otherwise  different  in  its  hierarchical  approach.  The  new  hierarchical  reconstruction  involves  a  MUSCL 
or  a  second  order  ENO  reconstruction  in  each  stage  of  a  multi-layer  reconstruction  process  without 
characteristic  decomposition.  It  is  compact,  easy  to  implement  over  arbitrary  meshes  and  retains  the 
overall  pre-processed  order  of  accuracy  while  effectively  removes  spurious  oscillations  around  shocks. 
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1.  Introduction 

The  first-order  Godunov  and  Lax-Friedrichs  (LxF)  schemes  are,  respectively,  the  forerunners  for  the 
large  classes  of  upwind  and  central  high-resolution  schemes  for  nonlinear  conservation  laws  and  related 
equations.  The  Godunov  scheme  captures  shock  waves  monotonically  in  narrow  transition  layers.  It  is 
based  on  evolving  a  piece- wise  cell  average  representation  of  the  solution,  by  evaluating  the  fluxes  at  the 
boundaries  of  each  cell  which  are  obtained  from  the  solution  of  (approximate)  Rienrann  problems  along 
the  boundary  interfaces.  Various  higher-order  generalizations  of  Godunov  scheme  have  been  developed 
since  the  mid  70s.  They  employ  higher-order  piece-wise  polynomials  which  are  reconstructed  from  the 
evolving  cell  averages  ‘in  the  direction  of  smoothness’.  We  mention  here  the  notable  examples  of  the 
high-resolution  upwind  FCT,  MUSCL,  TVD,  PPM,  ENO,  and  WENO  schemes  [6,  42,  17,  14,  18,  29] 
and  this  list  is  far  from  being  complete.  The  use  of  intricate  Rienrann  solvers  can  be  avoided  at  the 
expense  of  using  the  more  diffusive  LxF  scheme.  The  excessive  numerical  dissipation  can  be  reduced 
significantly,  however,  when  higher  order  piece- wise  polynomial  reconstructions  are  used  in  conjunction 
with  the  staggered  formulation  of  the  LxF  scheme.  The  central  scheme  of  Nessyahu  and  Tadrnor  (NT) 
[33]  provides  such  a  second-order  generalization  of  the  staggered  LxF  scheme.  It  is  based  on  the 
same  piece- wise  linear  reconstructions  of  cell  averages  used  with  upwind  schemes,  yet  the  solution  of 
(approximate)  Rienrann  problems  is  avoided.  High  resolution  generalizations  of  the  NT  scheme  were 
developed  since  the  90s  as  the  class  of  central  schemes  in  e.g.  [36,  2,  20,  19,  30,  4,  22,  1,  23,  24,  27], 
and  here  too,  the  list  is  far  from  being  complete.  The  relaxation  scheme  of  Jin  and  Xin  [21]  provides 
another  approach  which  leads  to  a  staggered  central  stencil  for  solving  nonlinear  conservation  laws. 

Being  free  of  the  (eigen-structure  of  -)  the  underlying  Rienrann  problems,  central  schemes  provide 
black-box  type  methods  for  the  approximate  solution  of  nonlinear  hyperbolic  conservation  laws  and 
other  closely  related  equations  [3].  Essentially,  one  only  needs  to  supply  the  flux  functions.  But  the 
staggered  high-order  central  schemes  of  order,  say,  r  >  1,  still  share  one  disadvantage  with  the  original 
LxF  scheme,  namely,  the  amplitude  of  their  numerical  viscosity  of  order  0((Ax)r+1  /At).  It  excludes 
the  use  of  small  tinre-steps,  At,  which  are  too  small  relative  to  the  spatial  grid  size  Ax.  The  problem 
lies  with  the  space-time  control  volumes  which  are  staggered  ‘Ax/2-  away’  from  each  other.  (Similar 
difficulties  occur  with  the  two  dimensional  conservative  front  tracking  method  which  was  overcome 
by  Glinrm  et  al.  in  [15],  using  space-time  cells  instead  of  rezoning).  This  problem  was  addressed 
by  Kurganov  and  Tadrnor  who  introduced,  in  [24],  a  new  type  of  central  schemes  whose  numerical 
viscosity  is  independent  of  0(1/ At)).  This  was  achieved  by  using  variable  control  volumes  so  that 
cells  are  staggered  only  ‘(P(At)-away’  from  each  other.  It  allows  implementation  of  central  scheme 
with  arbitrarily  small  time-step,  and  in  particular,  it  yields  their  semi-discrete  formulation  which  can 
be  conveniently  integrated  by  ODE  solvers,  e.g.,  the  Strong  Stability  Preserving  (SSP)  Runge-Kutta 
methods  of  [38],  consult  e.g.,  [16].  Similar  advantages  of  a  semi-discrete  formulation  can  be  achieved 
when  a  local  Lax-Friedrichs  building  block  is  used  over  non-staggered  meshes,  see  e.g.  Shu  and  Osher 
[38,  39]  and  Liu  and  Osher  [28].  The  upwind  and  central  schemes  mentioned  so  far,  share  one  thing 
in  common  —  they  evolve  one  piece  of  information  per  cell,  that  is,  the  cell  average.  Upwind  schemes 
use  Rienrann  solvers  while  central  schemes  use  simpler  quadrature  rules.  For  higher  accuracy,  they 
both  employ  piece-wise  polynomial  representation  of  the  solution  which  is  reconstructed  from  these 
cell  averages. 

In  [31],  Y.  Liu  introduced  an  alternative  technique  to  control  the  numerical  dissipation  of  central 
schemes.  The  main  idea  is  to  evolve  the  solution  over  overlapping  cells.  That  is,  two  sets  of  cell- 
averages  are  realized  over  interlacing  grids.  The  solution  is  then  represented  as  a  convex  combination 
-  an  ‘ O ( A f) -weight ed ’  combination  of  these  overlapping  cell  averages.  The  resulting  scheme  has  a 
numerical  viscosity  which  is  independent  of  0(1  /At)  and  as  such,  it  admits  a  semi-discrete  formulation 
which  can  be  integrated  using  SSP  methods.  The  use  of  overlapping  cells,  however,  is  fundamentally 
different  in  that  it  evolves  two  independent  quantities  for  each  given  cell,  that  is,  the  two  overlapping 
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sub-cell  averages.  The  use  of  overlapping  cells  opens  many  new  possibilities.  For  example,  instead 
of  the  usual  reconstructions  such  as  MUSCL  and  (W)ENO,  overlapping  cells  offer  a  more  efficient 
approach  for  high-resolution:  by  adding  the  two  subcell  averages,  we  recover  the  evolution  of  a  full 
cell  average,  where  by  taking  their  difference,  we  independently  evolve  an  approximate  slope,  rather 
than  reconstructing  it  from  neighboring  averages.  This  makes  feasible  the  use  of  central  discontinuous 
Galerkin  (DG)  approach  over  overlapping  cell,  following  the  series  of  works  of  Cockburn  and  Shu 
[11,  12,  13].  Thus,  in  particular,  the  use  of  overlapping  cells  yields  the  versatility  of  finite  element 
(Galerkin)  methods  which  can  be  easily  formulated  on  general  unstructured  meshes  with  any  formal 
order,  since  no  reconstruction  is  involved.  In  this  paper,  we  develop  the  staggered  central  DG  method 
for  solving  hyperbolic  conservation  laws.  Numerical  tests  are  performed  up  to  third  order  accuracy 
on  uniform  staggered  meshes  in  one  and  two  dimensions.  Stability  analysis  and  error  estimates,  and 
extension  of  the  method  for  time  dependent  and  steady  state  convection-diffusion  equations,  constitute 
ongoing  work  and  will  be  reported  in  the  future.  Here,  one  does  not  reconstruct  a  piecewise-polynomial 
representation  of  the  solution,  but  rather  it  is  part  of  the  evolution  of  higher-moments.  Still,  a  nonlinear 
limiting  procedure  is  necessary  to  reduce  spurious  oscillations  for  high  order  methods.  We  introduce 
here  such  a  general  non-oscillatory  procedure,  so-called  hierarchical  reconstruction  interesting  for  its 
own  sake,  which  is  closely  related  to  the  moment  limiters  of  Biswas  et  al.  [5]  and  to  the  earlier  work 
of  Cockburn  and  Shu  [11].  Since  this  limiting  procedure  requires  only  linear  reconstructions  using 
information  from  adjacent  cells  without  characteristic  decomposition,  it  can  be  easily  implemented 
for  any  shapes  of  the  cells  and  hence  is  practical  also  for  unstructured  meshes  or  even  dynamically 
moving  meshes  (e.g.  Tang  and  Tang  [40]),  although  we  do  not  pursue  it  in  this  paper.  We  refer  the 
reader  to  [2]  and  the  reference  therein  for  a  systematic  study  of  central  schemes  on  unstructured  grids 
using  the  framework  of  discontinuous  finite  elements. 

This  paper  is  organized  as  follows.  In  Sec.  2,  we  briefly  describe  the  central  finite  volume  scheme  on 
overlapping  cells  as  the  background.  The  natural  extension  to  central  DG  scheme  on  overlapping  cells 
is  discussed  in  Sec.  3.  In  subsection  3.1  we  study  the  numerical  convergence  rate  for  a  number  of  linear 
and  nonlinear  equations  having  smooth  solutions.  In  Sec.  4,  we  introduce  a  general  non-oscillatory 
hierarchical  reconstruction  procedure  and  use  it  as  a  post-processing  limiter  for  the  central  DG  scheme 
on  overlapping  cells  to  control  spurious  oscillations  in  the  presence  of  shocks.  Numerical  results  testing 
the  accuracy  of  the  proposed  schemes  are  included  in  Sections  3  and  4.  Additional  numerical  results 
are  presented  in  Sec.  5.  Concluding  remarks  and  a  plan  for  future  work  are  included  in  Sec.  6. 


2.  Central  Schemes  on  Overlapping  Cells 
Consider  the  scalar  one  dimensional  conservation  law 


(2.1) 


du  df(u) 
dt  dx 


( x ,  t)  G  TZ  x  (0,  T). 


Set  {xi  :=  xq  +  iAx},  let  Q+i/2  :=  [xi,xi+ 1)  be  a  uniform  partition  of  1Z  and  let  {U™+1/2}  denote 
the  set  of  approximate  cell  averages  t//+1/2  ~  (1/Ax)  Jc^^u(x,tn)dx.  Similarly,  we  set  D*  := 

[xi-i/2i  xi+i/2)  as  the  dual  partition  and  let  {V™}  denote  the  corresponding  set  of  approximate  cell 
average  V™  «  (1/Ax)  jD  u(x,tn)dx.  Starting  with  these  two  piecewise-constant  approximation1, 


X  Ui+l/21Ci+1/2(x)  and  EKiaW 


1Here  and  below,  ln(x)  denotes  the  characteristic  function  of  Q 
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we  proceed  to  compute  our  approximate  solution  at  the  next  time  level,  t”+1  :=  t”  +  A t”.  To  this 
end,  we  reconstruct  two  higher-order  non-oscillatory  piecewise-polynomial  approximations, 


un{x)  =  Y^Ui+1/2(x)lci+1/2(x ) 

i 


and  V\x )  =  ^2  Vi(x) lDi(x) 

i 


with  breakpoints  at  Xi,i  =  0,  ±1,±2,  and  respectively,  at  Xi+i/2,i  =  0,  ±1,  ±2,  ■  •  • .  These 
piecewise-polynomials  should  be  conservative  in  the  sense  that  Jc-+1/2  Un(x)dx  =  A xU™+i/2  and 

Id,  v"(  x)dx  =  A xV'j  for  all  j’s.  There  are  large  libraries  for  such  conservative,  accurate  and  non- 
oscillatory  reconstructions;  we  refer,  for  example,  to  the  second-order  example  of  MUSCL  [41],  the 
third-order  example  of  [30],  the  well-known  class  of  high-order  (W)ENO  reconstructions  [18,  37],  etc. 
Following  Nessyahu  and  Tadnror  [33] ,  the  central  scheme  associated  with  these  piecewise-polynomials 
reads 


(2.2a)  VT1 

(2.2b)  U£- 1/2 


Un{x)dx  —  f(Un+Hxi+l/2))  -  f(Un+HXi_ i/2)) 


j_  r  , .  At 

Ax 

1  /'  At” 

—  /  Vn(x)dx~  — 
Ax  Jc v  '  Ax 


i+1/2 


f(Vn+L2(xi+i))~f(Vn+Hxi)) 


To  guarantee  second-order  accuracy,  the  right-hand-sides  of  (2.2a),  (2.2b)  require  the  approximate 
values  of  Un+^(xj+i/ 2)  ~  u(xj+1/2,  tn+ 5)  and  Vn+^(xj)  «  u(xj,  tn+ 5)  to  be  evaluated  at  the  midpoint 
t  +  At”/2.  Replacing  the  midpoint  rule  with  higher  order  quadratures,  yields  higher  order  accuracy, 
e.g.,  [30,  4], 

The  central  Nessyahu-Tadnror  (NT)  scheme  (2.2)  and  its  higher-order  generalizations  provide  ef¬ 
fective  high-resolution  “black-box”  solvers  to  a  wide  variety  of  nonlinear  conservation  laws.  When  At 
is  very  small,  however,  e.g.,  with  At  =  0((A.x)2)  as  required  by  the  CFL  condition  for  convection- 
diffusion  equations  for  example,  the  numerical  dissipation  of  the  NT  schemes  becomes  excessively 
large.  The  excessive  dissipation  is  due  to  the  staggered  grids  where  at  each  time-step,  cell  averages 
are  shifted  A 2; / 2- away  from  each  other:  indeed,  at  the  extreme  of  f(u)  =  0,  the  central  scheme  (2.2) 
is  reduced  to  re-averaging  at  every  time  step.  To  address  this  difficulty,  Kurganov  and  Tadnror,  [24] , 
suggested  to  remove  this  excessive  dissipation  by  using  staggered  grids  which  are  shifted  only  O(At)- 
away  from  each  other.  This  amounts  to  using  control  volumes  of  width  O(At)  so  that  the  resulting 
schemes  admits  semi-discrete  limit  as  At  — »  0,  the  so  called  “central- upwind”  schemes  introduced  in 
[24]  and  further  generalized  in  [23].  Liu  [31]  introduced  another  modification  of  the  NT  scheme  which 
removes  its  0(1/ At)  dependency  of  numerical  dissipation.  In  this  approach,  one  takes  advantage  of 
the  redundant  representation  of  the  solution  over  overlapping  cells,  K™  and  U™+i/2-  The  idea  is 
to  use  a  0(At)-dependent  weighted  average  of  U™+\li  and  V"”.  In  fact  the  difference  between  them 
is  the  local  dissipation  error.  To  simplify  our  discussion,  we  momentarily  give  up  on  second-order 
accuracy  in  time,  setting  Un+  2  =  Un  and  Vn+i  =  Vn  in  (2.2a)  and  (2.2b).  The  resulting  first-order 
forward-Euler  formulation  of  the  new  central  scheme,  consult  figure  1,  reads 


(2.3a) 


F/+1 


(2.3b)  t/r+i/2 


Un(x)dx )  +  (1  -9)Vl 


At” 

Ax 


f(Un(xi+1/2))  -  f(Un(Xi_1/2)) 


Vn(x)dx^J 


+  (l-0)CC+l/2 


At” 

Ax 


f(Vn(xi+ 1))  -  f(Vn(Xi)) 


Here  0  :=  At”/Ar”  where  At”  is  an  upper  bound  for  the  time  step,  dictated  by  the  CFL  condition. 
Note  that  when  9  =  1,  (2.3a),  (2.3b)  is  reduced  to  the  first-order,  forward-Euler-based  version  of  the 
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A  B  C 

Figure  1:  (A)  NT  scheme;  (B)  ID  overlapping  cells;  (C)  overlapping  cells  create  self  similarity  for  the 
grid  over  time  and  allow  a  convex  combination  of  the  overlapping  cell  averages  to  control  the  numerical 
dissipation. 


NT  scheme  (2.2a),  (2.2b).  Moreover,  writing 


i-k  L  u"(x)d *) + o  -  ■ m  - v" + ^  (s  L  u"{x)dx  - 


and  recalling  that  At"  =  0( Ax)  due  to  the  CFL  restriction,  it  follows  that  the  local  dissipative  error 
now  has  a  pre-factor  of  order,  At",  and  hence  the  cumulative  error  will  be  independent  of  O(At).  The 
reduced  dissipation  allows  us  to  pass  to  a  semi-discrete  formulation:  subtracting  V™  and  t/”+1/2  from 
both  sides,  multiplying  by  and  then  passing  to  the  limit  as  At"  ->0we  end  up  with 


(2.4a)  jV,(n 


A  Tr 


L  Un{x)dx  ~ "  h  [/(C/n(x*+V2))  -  f(Un(Xi_  1/2))]  ■ 


(2.4b^C/m/2(t")  = 


Vn{x)dx-Uni+ 1/2 


i+1/2 


Ax 


[/(V”(ii+1))-/(V"(n))]. 


The  spatial  accuracy  of  the  semi-discrete  central  scheme  (2.4)  is  dictated  by  the  order  the  reconstruc¬ 
tion  Un(x )  and  Vn(x).  The  SSP  Runge-Kutta  methods  yield  the  matching  high-order  discretization 
in  time. 

We  conclude  this  section  by  quoting  [31]  regarding  the  non-oscillatory  behavior  of  the  central  scheme 
(2.4),  which  is  quantified  here  in  terms  of  the  Total- Variation-Diminishing  (TVD)  property,  e.g.,  [17]. 


Theorem  1.  Consider  the  central  schemes  (2.2)  and  (2.3)  which  are  set  with  the  same  initial  values 
Vi  and  Ui+l/2  at  t  =  tn.  If  the  NT  scheme  is  TVD,  then  so  is  the  central  scheme  (2.3). 

There  are  two  reconstruction  procedures  for  overlapping  cells:  one  is  the  standard  procedure  to 
reconstruct  the  two  classes  of  cell  averages  {Vi  :  i  =  0,  ±1,  ±2,  •  •  •}  and  {Ui+i/2  :  i  =  0,  ±1,  ±2,  •  •  •}; 
the  other  couples  these  two  classes  for  reconstruction  of  the  final  representation  of  the  solution.  Thus, 
this  approach  is  redundant.  At  the  same  time,  numerical  examples  in  [31]  have  shown  that  by  coupling 
the  reconstructions,  redundancy  does  provides  improved  resolution  when  compared  with  the  one-cell 
average  evolution  approach  of  Godunov-type  schemes. 


3.  A  Central  Discontinuous  Galerkin  Method  on  Overlapping  Cells  for 

Conservation  Laws 


Following  the  general  strategy  of  the  discontinuous  Galerkin  (DG)  methods,  see  e.g.  Lesaint 
and  Raviart  ([26]),  Cockburn  ([8]),  Cockburn  and  Shu  ([11,  13])  etc.,  the  central-type  discontinuous 
Galerkin  (DG)  method  on  overlapping  cells  can  be  derived  [32].  Consider  the  system  of  conservation 
laws 


duk 


+  Vx  •  ffc(u)  =  0, 


(x,  t)  S  IZd  x  (0,  T),  k  =  1, . . .  ,  m, 


(3.1) 


dt 
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\  / 


\ 

1 - 

x 

Figure  2:  2D  overlapping  cells  by  collapsing  the  staggered  dual  cells  on  two  adjacent  time  levels  to 
one  time  level. 

where  u  =  (it/, . . .  ,um)T.  For  simplicity,  assume  a  uniform  staggered  rectangular  mesh  depicted  in 
figure  2  for  the  2D  case,  and  we  note  that  the  similar  formulation  for  irregular  staggered  meshes,  e.g., 
the  Voronoi  mesh  consisting  of  a  triangular  mesh  and  its  dual. 

Let  {Cj+1/2},  I  =  •  •  •  ,  id)  be  a  partition  of  Rd  into  uniform  square  cells  depicted  by  solid 

lines  in  figure  2  and  tagged  by  their  cell  centroids  at  the  half  integers,  x/+1/2  :=  (I  +  1/2)  Ax.  Let  A4 
denote  the  set  of  piece- wise  polynomials  of  degree  r  on  the  cells  {CI+ 1/2};  no  continuity  is  assumed 
across  cell  boundaries.  Let  {Dj}  be  the  dual  mesh  which  consists  of  a  Ax/2-  shift  of  the  Cj+  j  /2 ’s 
depicted  by  dash  lines  in  figure  2.  Let  x/  be  the  cell  centroid  of  the  cell  Dj  and  let  AT  denote  the  set 
of  piece-wise  polynomials  of  degree  r  over  the  cells  {Dj};  again,  no  continuity  assumed  across  the  cell 
boundary.  The  weak  formulation  of  (3.1)  over  these  cells  reads 

d  f 

(3.2a)  —  /  uk(pdx 

M  JCI+ 1/2 

d  f 

(3.2b)  —  /  rtfci/’dx 

JD j 

where  n  is  the  unit  outer  normal  of  the  corresponding  cell  and  tp  and  ip  are  test  functions.  As  in  the 
one-dinrensional  setup,  we  let 

Un(x)  =  G  M  and  V”(x)  =  G  M 

1+ 1/2  i 

denote  two  representations  of  the  numerical  solution,  approximating  u(-,  tn)  over  the  two  overlapping 
grids,  {Cj-+1/2}  and  {Dj}.  Observe  that  each  of  the  two  vector  functions,  Un  with  smooth  pieces 
supported  on  the  C/+1/2’s  and  V”  with  smooth  pieces  supported  on  the  -D/’s,  consist  of  nn  components 

U"(x)  =  (UT(x), . . .  ,  t/;/(x))T  and  Vn(x)  =  (Vf  (x), . . .  ,  V£(x))T. 

Given  these  conservative,  accurate  and  non-oscillatory  approximations  at  tn  we  proceed  to  compute 
the  approximate  solution  at  the  next  time  level,  tn+1  =  tn  +  A tn.  To  this  end,  the  exact  solution  of 
the  right  hand  side  of  (3.2a),  is  replaced  by  Vn(x)  =  (V”, . . .  ,  V)”)T;  similarly,  for  the  right-hand-side 
of  (3.2b)  we  use  the  approximate  solution  Un(x)  =  (Df, . . .  ,  D,”)T.  Further,  time  derivatives  on  the 
left  are  replaced  by  Forward  Euler  time  differencing  where  we  use  the  same  type  of  0-weighting  of  the 


'  C/+ 1/2 


f k  ■  Vx(/dx  - 


/  (ffc  •  n )4>ds,  \/(j)  6  M,  k  =  l,---,m, 

J  dCI+i/2 

/  ffc  •  Vx^’dx  —  /  (ffc  •  n )ij}ds,  Vi/  e  AT,  k  =  l,---,m, 

J  Dj  JdD j 
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Un’s  and  the  Vn’s  as  in  (2.3a),  (2.3b).  In  the  resulting  central  discontinuous  Galerkin  method  one 
seeks  piece- wise  polynomials  {U”+*,2}  £  M  and  {Vj+1}  £  M  such  that  for  all  F s, 


t/"+1(x)0(x)dx  =  6 


/ C , 


7+1/2 


lc 


K*(x)<?Hx)dx  +  (1  -  9) 


7+1/2 


1C, 


7+1/2 


(3.3a) 


+Af1 


>c, 


ffc(Vn(x))  •  Vx0dx 


7+1/2 


yfen+1(x)V’(x)dx 


(3.3b) 


—Atn 


7+1/2 


(ffc(Vn(x))  •  n)0(x)ds, 


V(j>  £  A4,  k  =  1,  •  •  •  ,  m, 


of  Uk  (x)^(x)dx  +  (1  -  61)  [  Vk(x)i>(x)dx 
Jd !  Jdj 

+Atn  [  ffc(Un(x))  •  Vx^(x)dx 

JD! 


—At'1  f  (ffc(Un(x))  •  n)V’(x)ds,  Vi/’  £  A f,  k  =  1,  •  •  •  ,  m. 

J  dDi 

Here  0  =  Atn/Arn  <  1,  is  the  maximum  time  step  size  determined  by  the  CFL  restriction, 
Afn  =  tn+1  —  tn  is  the  current  time  step  size.  The  resulting  central  DG  is  the  two-dimensional 
analogue  of  the  one-dimensional  central  scheme  (2.4).  And  as  in  the  one-dinrensional  case,  serni- 
discrete  version  of  (3.3)  can  be  obtained;  higher  order  Runge-Kutta  time  discretization  can  be  used 
to  match  the  high-order  accuracy  of  the  spatial  reconstructions.  We  conclude  with  the  senri-discrete 
central  DG  approximation  of  (3.1)  such  that  for  all  admissible  test  functions  0  and  and  all  F s, 


d 

dt 


f  Uk(t>dx \tn  =  — f  (14n(x)  -  E/£  (x))0(x)dx 

J  Cj_ |_i/2  ^  Ctj-1  /o 


Cl+1/2 


(3.4a) 


+  /  ffc(Vn(x))  •  Vx0dx 
J  C7+1/2 


lac 


(ffc(Vn(x))  •  n)0(x)ds,  V0  £  M,  k  =  1,  •  •  •  ,  m, 


1+1/2 


—  /  Vk-i/jdxltn 

dt  JDt 


{Uk  (x)  -  ^T(x))^(x)dx 


1  Dr 


(3.4b) 


-I-  /  ffc(Un(x))  •  Vx^(x)dx 

JDj 


'8D ! 


(ffc(Un(x))  •  n)'0(x)ds,  \/ip  £  J\T,  k  =  1,  •  •  •  ,  m. 


3.1.  Numerical  Errors  For  Smooth  Solutions.  In  this  subsection  we  study  the  convergence  rate 
for  a  number  of  equations  having  smooth  solutions.  Let  us  start  with  the  following  linear  transport 
equation  with  periodic  boundary  conditions 

ut  +  ux  =  0,  (x,t)  £  (0,2)  x  (0,2) 


(3.5) 


u(x,  0)  =  1  +  sin(7r.T) 


x  £  (0,2) 
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Table  1.  P1 
(3.5). 


Ax 

1/20 

1/40 

1/80 

1/160 

1/320 

L\  error 

8.91E-3 

2.25E-3 

5.66E-4 

1.42E-4 

3.54E-5 

order 

- 

1.99 

1.99 

1.99 

Loo  error 

5.92E-3 

1.55E-3 

3.96E-4 

1.00E-4 

2.51E-5 

order 

- 

1.93 

1.97 

1.99 

1  version  of  t’ 

he  central 

DG  scheme  (3.4)  for  the  linear  convect 

Ax 

1/2 

1/4 

1/8 

1/16 

1/32 

L\  error 

6.69E-2 

3.29E-2 

5.04E-3 

1.66E-3 

3.88E-4 

order 

- 

1.02 

2.70 

1.60 

2.10 

Lqo  error 

3.85E-2 

2.05E-2 

7.69E-3 

1.19E-3 

2.75E-4 

order 

- 

0.91 

1.41 

2.69 

2.11 

Table  2.  P 1  version  of  the  central  DG  scheme  (3.4)  for  the  2D  Burgers’  equation 


(3.5). 


Ax 

1/20 

1/40 

1/80 

1/160 

1/320 

L\  error 

6.50E-5 

8.12E-6 

1.02E-6 

1.27E-7 

1.59E-8 

order 

- 

3.00 

2.99 

3.01 

Lqo  error 

4.68E-5 

5.90E-6 

7.40E-7 

9.27E-8 

1.16E-8 

order 

- 

2.99 

3.00 

3.00 

2  version  of  t’ 

he  central 

DG  scheme  (3.4)  for  the  linear  convect 

The  test  results  at  T  =  2  for  the  P 1  (piece- wise  linear)  version  of  the  central  DG  scheme  on 
overlapping  cells  (3.4)  are  listed  in  Table  1,  with  second  order  Runge-Kutta  time  discretization.  The 
CFL  factor  is  0.4  for  choosing  At  and  the  actual  time  step  size  At  is  chosen  with  8  =  0.9. 

It  can  be  seen  that  the  expected  second  order  accuracy  is  achieved.  We  also  conduct  a  convergence 
test  for  the  P 1  version  of  the  scheme  (3.4)  on  a  2D  problem  [10]  which  is  the  Burgers’  equation  with 
periodic  initial  data:  ut  +  {^u2)x  +  {\u2)y  =  0  on  [—1, 1]  x  [—1, 1],  u(x ,  y,  0)  =  \  +  ^  sin(7r(x  +  y)).  The 
numerical  solutions  are  computed  at  the  final  time  T  =  0.1  when  the  exact  solution  is  still  smooth. 
The  CFL  factor  is  0.4  for  choosing  Ar  and  the  actual  time  step  size  At  is  chosen  with  8  =  0.9.  The 
errors  are  shown  in  Table  2.  Again  we  observe  the  expected  second  order  accuracy.  The  results  for 
the  P 2  (piece-wise  quadratic)  version  of  the  scheme  (3.4)  for  the  linear  convection  equation  (3.5)  are 
listed  in  Table  3,  with  a  third  order  TVD  Runge-Kutta  time  discretization  [38].  We  can  see  that  the 
expected  third  order  accuracy  is  achieved. 

We  also  test  the  scheme  for  the  Burgers’  equation  ut  +  (4 u2)x  =  0,  u(x,  0)  =  4  +  4sin(7rx).  The 
errors  are  shown  in  Table  4  at  the  final  time  T  =  0.1  when  the  solution  is  still  smooth.  Further  we  test 
the  scheme  for  the  2D  Burgers’  equation  ut  +  (\u2)x  +  (\u2)y  =  0,  u(x,  0)  =  4  +  4  sin(7r(x  +  y)).  The 
errors  are  shown  in  Table  5  at  the  final  time  T  =  0.1.  The  solution  of  the  2D  Burgers’  equation  may 
contain  linear  waves,  hence  we  also  test  the  scheme  on  another  2D  equation  ut  +  (4 u2)x  +  {\uA)y  =  0, 
u(x,  0)  =  q  +  ^  sin(7r(x+y)).  The  accuracy  of  the  numerical  solution  is  shown  at  T  =  0.1  in  Table  6.  It 
seems  that  for  all  these  cases  the  expected  third  order  accuracy  is  achieved,  at  least  for  the  L\  errors. 
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Ax 

1/10 

1/20 

1/40 

1/80 

1/160 

L\  error 

2.72E-5 

3.41E-6 

4.29E-7 

5.37E-8 

6.78E-9 

order 

- 

3.00 

2.99 

3.00 

2.99 

Lqo  error 

4.00E-5 

7.06E-6 

8.27E-7 

1.04E-7 

1.31E-8 

order 

- 

2.50 

3.09 

2.99 

2.99 

Table  4.  P2  version  of  the  central  DG  scheme  (3.4)  for  the  ID  Burgers’  equation. 


Ax 

1/4 

1/8 

1/16 

1/32 

1/64 

Li  error 

8.33E-3 

9.58E-4 

1.36E-4 

1.65E-5 

2.14E-6 

order 

- 

3.12 

2.82 

3.04 

2.95 

Lqo  error 

4.56E-3 

8.20E-4 

1.48E-4 

1.95E-5 

2.58E-6 

order 

- 

2.48 

2.47 

2.92 

2.92 

Table  5.  P 2  version  of  the  central  DG  scheme  (3.4)  for  the  2D  Burgers’  equation. 


Ax 

1/4 

1/8 

1/16 

1/32 

1/64 

L\  error 

5.35E-3 

5.75E-4 

6.80E-5 

7.81E-6 

9.77E-7 

order 

- 

3.22 

3.08 

3.12 

3.00 

Lqo  error 

2.57E-3 

3.16E-4 

8.00E-5 

1.10E-5 

1.53E-6 

order 

- 

3.02 

1.98 

2.86 

2.85 

Table  6.  P 2  version  of  the  central  DG  scheme  (3.4)  for  the  2D  nonlinear  equation. 

4.  A  General  Non-Oscillatory  Hierarchical  Reconstruction  Procedure 

Compared  to  finite  volume  schemes  which  evolve  only  cell  averages  over  time,  discontinuous  Galerkin 
methods  compute  and  evolve  a  high  order  polynomial  in  each  cell.  How  to  take  advantage  of  the  extra 
information  provided  by  the  DG  method  in  each  cell  and  use  it  in  the  limiting  process  where  the 
solution  is  non-smooth  is  a  challenge.  The  first  idea  is  given  by  Cockburn  and  Shu  [11]  for  the  DG 
method  which  limits  the  variation  between  a  cell  edge  value  and  its  cell  average  by  the  differences 
between  the  cell  averages  of  the  current  and  neighboring  cells.  The  higher  Legendre  moments  are 
truncated  in  a  cell  if  non-smoothness  is  detected.  This  process  is  shown  to  be  total  variation  bounded 
in  the  means.  A  generalization  is  introduced  in  Biswas  et  al.  [5],  which  detects  the  non-snroothness 
in  higher  degree  moments  and  applies  the  limiting  when  necessary  from  higher  to  lower  moments. 
In  Qiu  and  Shu  [35,  34],  a  high  order  WENO  reconstruction  is  used  as  a  limiter  for  the  so-called 
“trouble”  cells,  where  the  polynomial  defined  at  Gaussian  points  is  reconstructed  from  a  WENO 
procedure  and  is  projected  back  to  the  finite  element  space  to  replace  the  one  computed  by  the  DG 
method.  In  [34],  the  Hermite  WENO  reconstruction  takes  not  only  cell  averages  of  a  function,  but 
also  cell  averages  of  its  first  order  derivatives  in  order  to  obtain  a  compact  reconstruction.  A  similar 
strategy  is  used  in  our  non-oscillatory  hierarchical  reconstruction,  where  cell  averages  of  various  orders 
of  derivatives  of  a  function  are  to  be  calculated  and  used  in  the  reconstruction  of  linear  polynomials 
in  each  stage.  Our  limiting  procedure  is  closely  related  to  [5] .  Our  departure  from  [5]  begins  with  a 
different  point  of  view,  where  the  approximation  in  each  cell  is  viewed  as  a  high  degree  polynomial, 
instead  of  the  combination  of  orthogonal  Legendre  polynomials  advocated  in  [5] .  Instead  of  a  limiting 
procedure  which  is  trying  to  set  an  acceptable  range  for  the  coefficient  of  the  Legendre  moments 
(by  using  the  coefficients  of  lower-degree  moments),  we  reconstruct  the  complete  set  of  coefficients 
of  the  m-degree  polynomial  terms,  using  a  non-oscillatory  conservative  reconstruction  which  involves 
previous  reconstructed  terms  of  degrees  above  m.  The  resulting,  so-called  hierarchical  reconstruction 
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algorithm  is  easy  to  implement  in  a  multi-dimensional  setting  and  there  is  no  need  to  transform  an 
irregular  mesh  cell  into  a  rectangular  one  or  use  a  dimension-by-dimension  extension  of  a  ID  limiter. 
It  is  essentially  independent  of  the  shapes  of  mesh  cells  and  is  compact  because  of  the  conservative 
non-oscillatory  linear  reconstruction  (such  as  the  MUSCL  or  second  order  ENO  reconstruction)  used 
in  each  stage.  We  now  give  the  details  of  this  reconstruction  procedure.  For  simplicity  we  only  discuss 
the  scalar  case.  For  systems  a  component-by-component  extension  is  applied  without  characteristic 
decomposition. 

From  scheme  (3.4)  with  the  SSP  Runge-Kutta  methods,  we  obtain  numerical  solutions  Un(x)  and 
Un(x)  at  the  time  tn.  To  simplify  the  notation  we  drop  the  superscript  n  and  write 

u(*)  =  Y  Ul+ V2(x  “  x/+i/2)1C7+1/2(x)  G  M  and  V(x)  =  Y  VHX  ~  x/)1w(x)  e  ^ 

1+ 1/2  / 


recalling  that  x/+1/2  and  x/  are  centroids  of  cell  Cj+ \/2  and  Dj  respectively;  Uj+ i/2(x  —  x/+1/2)  and 
V/(x  —  x/)  are  the  polynomials  (of  degree  r)  in  cells  Cj+X  i2  and  Dj  respectively2.  The  task  is  to 
reconstruct  a  ’limited’  version  of  these  polynomials,  retaining  high-resolution  and  removing  spurious 
oscillations.  In  the  following,  we  discuss  a  hierarchical  reconstruction  procedure  to  recompute  the 
polynomial  UI+ i/2(x  —  x/_ (_1/2)  by  using  polynomials  in  cells  adjacent  to  cell  Cj+i/2-  For  convenience 
these  adjacent  cells  are  renamed  as  the  set  {Cj}  (which  contain  cell  C/+1/2,  Dj  etc),  and  the  poly¬ 
nomials  (of  degree  r)  supported  on  them  are  thus  renamed  as  {I/j(x  —  xj)}  respectively,  where  xj  is 
the  cell  centroid  of  cell  Cj.  We  write  Uj+ i/2(x  —  x/+1/2)  in  terms  of  its  Taylor  expansion, 

UI+l/2(x  -  X/+1/2)  =  Y  Y  ^!^S/2(°)(x-xm/2)m, 

in= 0  |m|  =m 


where  _ 

m!  /+1/2 


(0)  are  the  coefficients  which  participate  in  its  typical  m-degree  terms, 


Y  lml  =  0, .  . .  ,  r. 

|m|  =m 

The  following  hierarchical  reconstruction  describes  a  procedure  to  compute  the  new  coefficients 

~[^/+i/2(0)’  m  =  r,r-  1,...,0 

in  Uj_ |_i/2(x  —  xj+1/2),  iterating  from  the  highest  to  the  lowest  degree  terms. 


4.1.  An  Example  for  Piece-wise  Quadratic  Finite  Element  Space  in  ID.  Suppose  Uj(x  — 
Xj)  =  Uj( 0)  +  Uj(0)(x  —  Xj)  +  \Uj{ 0)(x  —  Xj)2,  j  =  1,2,3,  are  given  at  cells  C\,  C2  and  C3  respectively 
(see  figure  3  left),  where  Xj  is  the  center  of  cell  Cj.  These  polynomials  could  be  oscillatory  if  located 
near  a  discontinuity  of  the  weak  solution.  The  following  algorithm  computes  a  new  coefficient  for  each 
coefficient  in  the  polynomial  defined  on  cell  C2,  in  order  to  reduce  the  oscillation  while  keeping  the 
accuracy  (in  smooth  area)  and  resolution. 

First  take  the  first  derivative  for  them  to  obtain  Lj(x  —  Xj)  =  U'j( 0)  +  U"{ 0)(x  —  Xj).  Then  calculate 
the  cell  average  of  Lj(x  —  Xj)  on  cell  Cj  to  obtain  Lj  =  Uj( 0),  j  =  1,2,3.  With  the  three  cell 
averages  one  can  apply  a  MUSCL  or  second  order  ENO  procedure  to  reconstruct  a  non-oscillatory 
linear  polynomial  in  cell  C2.  For  the  MUSCL  reconstruction  with  the  minrnod  limiter,  we  obtain  a 
new  coefficient  for  U2  (0),  denoted  by 


^'(0) 


minrnod 


L2  —  L\ 

x2  -  X\ 


L3  —  L2 
x3  ~  x2 


2These  polynomials  could  be  oscillatory.  There  could  be  other  methods  to  compute  these  polynomials  such  as  a  finite 
volume  reconstruction  from  cell  averages. 
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Figure  3:  Left:  ID  non-oscillatory  hierarchical  reconstruction  for  cell  2  involves  only  overlapping  cells 
1,2  and  3.  Right:  2D  non-oscillatory  hierarchical  reconstruction  for  cell  3  involves  only  overlapping 
cells  1, 2,  3, 4  and  5. 


where 


{min{ci,c2,  •  •  -  ,  cm},  ifci,c2  5***5  Cm  >  o, 
max{ci,c2,-  ■  •  ,cm},  ifci,c2 

5***5  Cm  <  o, 

0,  otherwise. 


For  the  second  order  ENO  reconstruction, 


U2  (0)  =  minmod2 


L2  —  Li  L%  —  L2 

X2~X\'  X3  -  X2 


where 


minmod2{ci,  c2,  •  •  •  ,  cm}  =  Cj ,  if  \c3\  =  min{|ci|,  |c2|,  •  •  •  ,  |cm|}. 


In  order  to  find  a  new  coefficient  U2(0)  for  U2( 0),  we  first  compute  the  cell  average  of  Uj(x  —  Xj)  on  cell 
Cj  to  obtain  Uj,  and  then  compute  the  cell  average  of  the  polynomial  i?2(x  — x2)  =  \U2 (0)(x  —  x2)2  on 
cell  Cj  to  obtain  Rj,  j  =  1, 2,  3.  Redefine  Lj  =  Uj  —  Rj,  j  =  1,2,3.  A  similar  MUSCL  or  second  order 
ENO  procedure  can  be  applied  to  the  three  cell  averages  L\ ,  L2  and  L3  to  find  U'2{ 0).  For  example, 
by  using  the  MUSCL  reconstruction,  we  obtain 


U2(0)  =  minrnod 


f  L2  —  Li  L3  —  L2 

l  X2  -  Xi  ’  X'3  -  X2 


Finally,  let  U2{ 0)  =  L2. 

The  convergence  test  results  with  Algorithm  1  for  the  same  problem  as  that  in  Table  4  can  be  found 
in  Tables  7  and  8.  We  observe  that  the  order  of  accuracy  is  maintained,  although  (as  expected  for  any 
limiter)  the  magnitude  of  the  error  is  increased  for  the  same  mesh. 


4.2.  Hierarchical  Reconstruction  —  General  Description.  In  the  following,  we  discuss  a  hierar¬ 
chical  reconstruction  procedure  to  recompute  the  polynomial  U/+]y2(x  —  x/+1/2)  by  using  polynomials 
in  cells  adjacent  to  cell  Cj+i/2.  Recall  that  these  adjacent  cells  are  renamed  as  the  set  {Cj}  and  the 
polynomials  (of  degree  r)  supported  on  them  are  thus  renamed  as  {Lj(x  —  xj)}  respectively.  The 
following  hierarchical  reconstruction  describes  a  procedure  to  compute  the  new  coefficients 

~[^/+i/2(0)’  m  =  r,r-  1,...,0 

in  UI+lj2(x  —  xJ+1/2),  iterating  from  the  highest  to  the  lowest  degree  terms. 
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To  reconstruct  R/+i/2(0)>  we  first  compute  many  candidates  of  ^j™i/2(0)  (sometimes  still  denoted 
as  U%}/2{0)  with  specification),  and  we  then  let  the  new  coefficient  for  t/j™ w2(0)  be 

^7+1/2 (0)  =  candidates  of  0)  ), 

where  F  is  a  convex  limiter  of  its  arguments,  e.g.,  the  minrnod  function. 

In  order  to  find  these  candidates  of  f/j™ jy2(0),  |m|  =  m,  we  take  a  (m  —  l)-th  order  partial  derivative 
of  hr/+i/2(x  —  x/+1/2),  and  denote  it  by 

d'"  1C//+1/2(x  -  x/+1/2)  =  T/+1/2(x  -  x/+1/2)  +  i?7+i/2(x  -  x/+1/2), 

where  Lj+i/2  is  the  linear  part  and  Rj. (_i/2  contains  all  second  and  higher  degree  terms  of  dm~1Uj+ 1 /2(x— 
x7+i/2).  Clearly,  every  coefficient  in  the  first  degree  terms  of  Tj+i/2  is  in  the  set  {^/™i/2(0)  :  |m|  =  m}. 
And  for  every  m  subject  to  |m|  =  rn,  one  can  always  take  some  (m  —  l)-th  order  partial  derivatives 
of  C/j_|_i/2(x  —  x/+1/2)  so  that  u£}/2{0)  is  a  coefficient  in  the  first  degree  terms  of  LI+1/2.  Thus,  a 
‘candidate’  for  a  coefficient  in  the  first  degree  terms  of  LI+l/2  is  the  candidate  for  the  corresponding 

O). 

In  order  to  find  the  candidates  for  all  the  coefficients  in  the  first  degree  terms  of  LI+l  /2(x  —  x/+1/2), 
we  only  need  to  know  the  new  approximate  cell  averages  of  T/+]y2(x— xj+1 /2)  on  d+ 1  distinct  mesh  cells 
adjacent  to  cell  Cj+\  /2,  recalling  that  d  is  the  spatial  dimension.  Assume  Cj0,  Cj1 ,  •  •  •  ,  Cjd  €  {Cj} 
are  these  cells  and  Lj0,  Ljx ,  •  •  •  ,  Ljd  are  the  corresponding  new  approximate  cell  averages.  The  set  of 
these  d+  1  cells  with  given  new  approximate  cell  averages  is  called  a  stencil.  Let  a  linear  polynomial 
Lm/2(x  —  xm/2)  be  determined  by 

(4.1)  W7]  Jc  ^/+1/2(x  ”  xm/2)^x  =  LJp  /  =  0,1,- ••  ,d. 

Then  the  coefficients  in  the  first  degree  terms  of  L/+i/2(x  —  xj+1/2)  will  be  the  candidates  for  the 
corresponding  coefficients  of  Lj+1/2(x  —  x/+1/2).  Therefore,  a  stencil  located  near  cell  Cj+i/2  will 
determine  a  set  of  candidates  for  all  coefficients  in  the  first  degree  terms  of  L/+1/2(x  —  x/+1/2).  The 
key  is  to  determine  the  new  approximate  cell  averages  of  Lj+1/2(x  —  x/+1/2)  on  the  cells  of  {Cj}, 
which  is  outlined  by  the  following  algorithm. 


Algorithm  1.  Step  1.  Suppose  r  >  2.  For  m  =  r,r  —  1,  •  •  •  ,  2,  do  the  following: 

(a)  Take  a  (m  —  1  )-th  order  partial  derivative  for  each  of  {Uj(x  —  xj)}  to  obtain  polynomials 
{9™-i t/j(x  — xj)}  respectively.  In  particular,  denote  dm~1Uj+i/2('x  —  x/+1/2)  =  L/+1/2(x  —  x/+1/2)  + 
4?/+i/2(x  -  x7+i/2)>  where  L/+1/2(x  -  xJ+1/2)  is  the  linear  part  of  5m_1C//+1/2(x  -  x/+1/2)  and 
i?/+i/2(x  —  x/+1/2)  is  the  remainder. 

(b)  Calculate  the  cell  averages  of  {dm~1Uj('x  —  xj)}  on  cells  { Cj }  to  obtain  {d'n~1Uj}  respectively. 

(c)  Let  ^7+1/2(x— x/+1/2)  be  the  -R/+i/2(x— x/+1/2)  with  its  coefficients  replaced  by  the  corresponding 
new  coefficients.  Calculate  the  cell  averages  of  RI+ i/2(x  —  x/+1/2)  on  cells  {Cj}  to  obtain  {-Rj} 
respectively. 

(d)  Let  Lj  =  dm~1Uj  —  Rj  for  all  J. 

(e)  Form  stencils  out  of  the  new  approximate  cell  averages  {Lj}  by  using  a  non- os  dilatory  finite 
volume  MUSCL  or  second  order  ENO  strategy.  Each  stencil  will  determine  a  set  of  candidates  for 
the  coefficients  in  the  first  degree  terms  of  L/+i/2(x  —  x/+1/2),  which  are  also  candidates  for  the 

corresponding  £/j™y2(0)  ’s,  |m|  =  m. 
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(f)  Repeat  from  (a)  to  (e)  until  all  possible  combinations  of  the  (m  —  1  )-th  order  partial  derivatives 
are  taken.  Then  the  candidates  for  all  coefficients  in  the  m-th  degree  terms  of  UI+ i/2(x  —  x/+ 1/2) 

have  been  computed.  For  each  of  these  coefficients,  say  j^T^j™i/2(0),  |m|  =  rn,  let  the  new  coefficient 

U%},2(0)  =  F(  candidates  of  ^j™i/2(0)  )• 

Step  2.  In  order  to  find  the  new  coefficients  in  the  zero-th  and  first  degree  terms  of  UI+ 1  /2 (x  ~ 
x/_i_i/2),  we  perform  the  procedure  of  Step  1  (a)-(f)  with  m  =  1,  and  make  sure  that  the  new 
approximate  cell  average  Lj_ j_i /2  is  in  each  of  the  stencils,  which  ensures  that  the  cell  average  of 
Ui+ i/2(x  —  x/+1/2)  on  cell  C/+1/2  is  not  changed  with  new  coefficients.  The  new  coefficient  in  the 
zero-th  degree  term  of  UI+ i/2(x  —  xj+1/2)  is  Lj+1/2. 

Remarks.  1.  The  coefficients  of  the  polynomials  can  be  updated  after  Algorithm  1  has  been  applied 
to  all  mesh  cells,  or  at  the  m-th  stage  when  all  new  coefficients  for  those  in  the  m-th  degree  terms  of 
all  polynomials  have  been  computed  (in  this  case,  {d°Uj}  used  in  Step  2  should  be  the  cell  averages 
of  the  original  polynomials  to  ensure  that  they  are  invariant).  The  latter  case  is  supposed  to  be  more 
diffusive.  In  numerical  experiments  we  find  their  results  are  about  the  same.  All  numerical  results 
presented  in  this  paper  are  performed  with  the  former  implementation. 

2.  One  motivation  for  us  to  develop  this  hierarchical  reconstruction  is  that  the  limiting  for  the  DG 
scheme  on  non-staggered  meshes  is  different  from  that  for  scheme  (3.4).  Because  for  the  usual  DG 
scheme  the  time  evolution  of  the  cell  averages  is  completely  determined  by  the  fluxes.  However  in 
(3.4),  cell  interior  values  are  also  involved.  We  find  in  numerical  experiments  that  the  moment  limiter 
[5]  does  not  work  as  robustly  for  scheme  (3.4)  as  it  does  for  DG  scheme  on  non-staggered  meshes.  The 
proposed  hierarchical  reconstruction  process  is  quite  general  and  could  be  useful  for  conventional  DG 
or  even  finite  volume  schemes.  These  will  be  explored  in  the  future. 

3.  Scheme  (3.4)  with  Algorithm  1  and  with  piece-wise  linear  elements  is  identical  to  the  second 
order  central  scheme  on  overlapping  cells  [31]. 

4.  It  is  more  efficient  to  apply  the  hierarchical  reconstruction  process  only  at  places  where  it  is 
needed  by  using  non-snroothness  detectors  (see  e.g.  [35,  7]).  This  will  be  explored  in  the  future. 

The  most  important  point  is  that,  even  though  the  linear  reconstruction  used  in  Algorithm  1  is  only 
second  order  accurate,  the  approximation  order  of  accuracy  of  a  polynomial  in  a  cell  is  unaffected  by 
the  algorithm  and  we  have  the  following  theorem. 

Condition  1.  Let  {xj0 ,  x  ,  •  •  •  ,xjd}  be  the  d  +  1  cell  centroids  of  a  stencil.  Then  there  is  a  point 
among  them,  say  xj0,  such  that  the  matrix  A  =  ^[xjt  — xj0,  xj2  — xj0,  •  •  •  ,  xjd  —  xj0]  is  non  singular. 
Further,  there  is  a  constant  a  >  0  independent  of  Ax  such  that  ||A_1||  <  a. 

Theorem  2.  Suppose  {Uj(x  —  x j)}  in  Algorithm  1  approximate  a  Cr+l  function  u(x)  with  point- 
wise  error  0((Ax)r+1)  within  their  respective  cell  {Cj},  and  all  cells  in  { Cj }  are  contained  in  a  circle 
centered  at  xJ+1/2  with  radius  0( Ax).  Let  the  d- 1-1  cell  centroids  in  every  stencil  used  in  Algorithm  1 
satisfy  Condition  1.  Then  after  the  application  of  Algorithm  1,  the  polynomial  D/+1/2(x  —  xj+1/2),  i.e. 
^/+i/2(x  —  x/+i/2)  its  coefficients  replaced  by  the  corresponding  new  coefficients  also  approximates 

the  function  u(x)  with  point-wise  error  C((Ax)r+1)  within  cell  C/+1/2.  The  cell  average  of  Uj+ i/2(x  — 
x/+ 1/2)  on  ce/Z  C/+1/2  ts  the  same  as  that  of  Uj+ i/2(x  —  xJ+1/2). 

Proof.  From  the  assumption  we  know  that  the  coefficients  in  the  m-th  degree  terms  of  £7/+i/2(x  — 
x/+ 1/2),  0  <  m  <  r,  are  the  (r  —  m  +  l)-th  order  approximation  to  the  corresponding  coefficients  of 
the  Taylor  expansion  of  u(x)  at  xI+l  /2. 

Assume  that  when  starting  to  compute  new  coefficients  for  coefficients  of  the  m-th  degree  terms 
of  UI+1/2fix  —  x/+1/2),  1  <  m  <  r,  all  computed  new  coefficients  (if  there  is  any)  for  coefficients 
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Ax 

1/10 

1/20 

1/40 

1/80 

1/160 

L\  error 

4.24E-4 

5.33E-5 

6.71E-6 

8.44E-7 

1.07E-7 

order 

- 

2.99 

2.99 

2.99 

2.98 

Lqo  error 

5.13E-4 

6.20E-5 

7.38E-6 

1.29E-6 

2.66E-7 

order 

- 

3.05 

3.07 

2.52 

2.28 

Table  7.  P 2  version  of  the  central 


DG  scheme  (3.4)  with  the  hierarchical  reconstruc¬ 


tion  Algorithm  1  for  the  Burgers’  equation.  MUSCL  is  used  in  Algorithm  1. 


of  the  1-th  degree  terms  (m  <  l  <  r,  if  they  exist)  of  C7/+]y2(x  —  xj+1/2)  are  their  (r  —  l  +  l)-th 
order  approximations.  In  fact,  when  m  =  r,  there  is  no  new  coefficients  been  computed  at  Step  1 
(a).  However,  the  following  argument  will  show  that  the  new  coefficients  computed  at  Step  1  (f)  for 
coefficients  of  the  ?’-th  degree  terms  of  H/+1/2(x  —  ^-1+1/2)  are  their  first  order  approximations. 

Let  Lj+1/2(x  — x/+1/2)  =  co+ci  -(x  — x/+1/2)  in  Step  1  (a)  and  let  L(x)  =  co+ci  -(x  —  x/+1/2)  be  the 
corresponding  linear  part  in  the  Taylor  expansion  of  the  same  (  as  for  Uj  )  (m  —  l)-th  partial  derivative 
of  u(x)  at  xJ+1/2.  Therefore  Co  and  ci  approximate  Co  and  ci  to  the  order  of  0[( Ax)r~m+2)  and 
0(( Ax)r_rn+1)  respectively.  Also  from  the  above  assumptions  it  is  easy  to  see  that  Lj  =  dm~1Uj  —  Rj 
in  Step  1  (d)  approximates  the  cell  average  of  L(x)  on  cell  Cj  to  the  order  of  0( Axr~m+2),  for  all  cell 
Cj  adjacent  to  cell  Cr+l  /2. 

Reconstructing  L/+1/2(x  -  x/+1/2)  =  c0  +  Ci  •  (x  -  xm/2)  from  a  stencil  Cj0,Cji:---  ,  Cjd  G  { Cj } 
is  to  find  Co  and  ci  satisfying  the  following  equations  (see  (4.1)), 

(4-2)  Jf=r~\jc  (co  +  ci  •  (x  -  x/+1/2))dx  =  c0 +  ci  •  (xj,  -  x/+1/2) 

=  LJt  =  c0  +  ci  •  (xJ;  -  xm/2)  +  £>(( Ax)r~m+2), 

where  x^  is  the  cell  centroid  of  cell  Cjt .  I  =  0,  •  •  •  ,d.  The  solutions  are  candidates  for  co  and  ci 
respectively.  Subtracting  the  first  equation  (/  =  0)  from  the  rest  of  the  equations  in  (4.2)  we  can 
obtain 

AT(c1-c1)  =  0((AxY~m+1), 

where  A  =  [xjj  —  xj0,xj2  —  xj0,  •  •  •  , xjd  —  xj0].  From  Condition  1,  1 1 1 1  is  bounded  independent 
of  Ax.  We  conclude  that  the  candidate 

(4.3)  C!  =C!  +  0((Ax)r-m+1). 

Also  since  ||xj(  —  x/+1/2||  =  O(Ax),  l  =  0, 1,  ■  ■  ■  ,  d,  by  substituting  the  estimate  of  the  candidate  Ci 
back  into  one  of  the  equations  of  (4.2)  we  obtain  that  the  candidate 

(4.4)  c0  =  c0  +  O{(AxY~m+2). 

Since  the  function  F  used  in  Step  1  (f)  is  a  convex  combination  of  its  arguments,  it  does  not  change 
the  approximation  order  of  its  arguments.  Therefore  estimate  (4.3)  implies  that  the  new  coefficients 
for  coefficients  of  the  m-th  degree  terms  of  UI+ i/2(x  —  xj+1/2)  are  their  (r  —  m  +  l)-th  order  approx¬ 
imations.  Estimate  (4.3)  moves  the  induction  till  m  =  1  and  estimate  (4.4)  implies  that  in  Step  2, 
the  new  coefficient  for  the  coefficient  of  the  zero-th  degree  term  of  L/+1/2(x  —  x/+1/2)  is  its  0(Axr+1) 
approximation.  Step  2  clearly  insures  that  the  cell  average  of  t//+1/2(x  —  x/+1/2)  on  cell  C/+1/2  is 
unchanged  with  new  coefficients.  The  proof  is  complete. 
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Ax 

1/10 

1/20 

1/40 

1/80 

1/160 

L\  error 

4.51E-4 

5.36E-5 

6.85E-6 

8.54E-7 

1.08E-7 

order 

- 

3.07 

2.97 

3.00 

2.98 

Lqo  error 

5.24E-4 

6.17E-5 

1.03E-5 

1.81E-6 

3.27E-7 

order 

- 

3.09 

2.58 

2.51 

2.47 

Table  8.  P 2  version  of  the  central 


DG  scheme  (3.4)  with  the  hierarchical  reconstruc¬ 


tion  Algorithm  1  for  the  Burgers’  equation.  Second  order  ENO  is  used  in  Algorithm  1. 


Ax 

1/4 

1/8 

1/16 

1/32 

1/64 

Li  error 

8.00E-2 

1.24E-2 

1.58E-3 

1.92E-4 

2.40E-5 

order 

- 

2.69 

2.97 

3.04 

3.00 

Lqo  error 

4.90E-2 

9.85E-3 

1.68E-3 

2.01E-4 

2.68E-5 

order 

- 

2.31 

2.55 

3.06 

2.91 

Table  9.  P 2  version  of  the  centra’ 
tion  Algorithm  1  for  the  2D  Burgers’  equation, 
rithrn  1. 


DG  scheme  (3.4)  with  the  hierarchical  reconstruc- 
Second  order  ENO  is  used  in  Algo- 


4.3.  Implementation  for  Piece-wise  Quadratic  Finite  Element  Space  in  2D.  Suppose  on  cell 
Cj  (see  figure  3  right),  a  quadratic  polynomial  is  given  as 

Uj(x  -  Xj,y  —  yj)  =  Uj( 0,0)  +  dxUj( 0,0) (a;  -  xj )  +  dyUj(0,0)(y  -  yj)  + 

^dxxUj( 0,  0)(x  -  Xj)2  +  dxyUj( 0,  0)(x  -  xj)(y  -  yj)  +  ^dyyUj(0,  0 )(y  -  y5 )2, 

where  (xj,  yj)  is  the  cell  centroid  of  cell  Cj,  j  =  1,  2,  •  •  •  ,5. 

According  to  Step  1  of  Algorithm  1,  take  the  first  partial  derivative  with  respect  to  x  for  them 
to  obtain  Lj(x  -  Xj,y-  yj)  =  dxUj( 0,0)  +  dxxUj(0,0)(x  -  Xj)  +  dxyUj(0,  0)(y  -  yj),  j  =  1,2,  •  •  •  ,5. 
Calculate  the  cell  average  of  Lj(x  —  Xj,  y  —  yj)  on  cell  Cj  to  obtain  Lj  =  dxUj(0,  0),  j  =  1, 2,  •  •  •  ,  5 
(note  that  R3(x  —  X3,  y  —  y3)  =  0).  With  the  five  new  approximate  cell  averages  {Lj  :  j  =  1, 2,  •  •  •  ,  5}, 
one  can  apply  a  MUSCL  or  a  second  order  ENO  procedure  to  reconstruct  a  non-oscillatory  linear 
polynomial 

L3(x  -x3,y-  2/3)  =  dxU3(0,  0)  +  dxxU3(0,  0)(x  -  x3)  +  dxyU3( 0,  0 )(y  -  y3) 

in  cell  C3.  For  example,  one  can  form  four  stencils  {C3,C\,C2},  {C3,  C2,  C5},  {C3,C3,Ci}  and 
{C3,  C4,  Ci}.  For  the  first  stencil,  solve  the  following  equations  for  dxxU3( 0,  0)  and  dxyU3( 0,  0) 

T7T1  [  L3(x  -  x3,y  -  y3)dxdy  =  L3  +  dxxU3(0,  0)(xj  -  x3)  +  dxyU3(0,  0)(yj  -  y3) 
l°il  JCj 

—  Lj, 

j  =  1,2,  similarly  for  other  stencils.  We  obtain  two  sets  of  candidates  for  dxxU3( 0,  0)  and  dxyU3(0,  0) 
respectively.  By  taking  the  first  partial  derivative  with  respect  to  y  for  Uj(x—Xj,  y—yj),  j  =  1,2,---  ,5, 
we  similarly  obtain  a  set  of  candidates  for  dyyU3( 0,  0)  and  enlarge  the  set  of  candidates  for  dxyU3((),  0). 
Putting  all  candidates  for  dxxU3( 0,0)  into  the  arguments  of  the  minrnod  (or  minimi)  function,  we 
obtain  the  new  coefficient  dxxU3((),  0)  for  dxxU3(0, 0).  Applying  the  same  procedure  to  obtain  new 
coefficients  dxyU3( 0,0)  and  dyyU3( 0,0). 
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According  to  Step  2  of  Algorithm  1,  we  compute  the  cell  average  of  Uj(x  —  Xj,y  —  yj)  on  cell  Cj  to 
obtain  Uj,  j  =  1,  2,  •  •  •  ,  5;  and  compute  cell  averages  of  the  polynomial 

Rz{x  -x3,y-  2/3 )  =  -dxxU3{ 0,  0)(x  -  x3 )2  +  dxyU3(0,  0)(x  -  x3)(y  -  y3 )  +  ^dyyU3(0, 0 ){y  -  y3 )2 

on  cell  C\ ,  C‘2 ,  •  •  ■  ,  C3  to  obtain  R\,  R2,  ■  ■  ■  ,  R3  respectively.  Redefine  Lj  =  Uj  —  Rj,  j  =  1,  2,  •  •  •  ,5. 
The  same  MUSCL  or  second  order  ENO  procedure  as  described  previously  can  be  applied  to  the  five 
cell  averages  {Lj  :  j  =  1,  2,  •  •  •  ,  5}  to  obtain  new  coefficients  dxU3(0,  0)  and  dyU3(0.  0).  Finally  let  the 
new  coefficient  U3(0,  0)  =  L3. 

The  convergence  test  results  with  Algorithm  1  for  the  same  problem  as  that  in  Table  5  can  be  found 
in  Table  9.  We  again  observe  that  the  order  of  accuracy  is  maintained,  although  (as  expected  for  any 
limiter)  the  magnitude  of  the  error  is  increased  for  the  same  mesh. 

5.  Additional  Numerical  Examples 

Scheme  (3.4)  with  the  piece-wise  r-th  degree  polynomial  space  is  referred  to  as  CO-DG-(r  +  1). 
When  the  hierarchical  reconstruction  Algorithm  1  is  applied,  it  is  referred  to  as  CO-DG-hrl-(r  +  1). 
To  specify  whether  a  linear  MUSCL  (with  the  minrnod  limiter)  or  ENO  (with  the  minmod2  limiter) 
reconstruction  is  used  in  Algorithm  1,  we  refer  it  as  CO-DG-hrlm-(r  +  1)  or  CO-DG-hrle-(?’  +  1) 
respectively. 

The  corresponding  (up  to  third  order)  TVD  Runge-Kutta  time  discretization  methods  [38]  are 
applied  to  the  above  schemes.  Only  the  solution  in  one  class  of  the  overlapping  cells  is  shown  in  the 
graphs  throughout  this  section.  For  systems  of  equations,  the  component-wise  extensions  of  the  scalar 
schemes  (without  characteristic  decomposition)  have  been  used  in  all  the  computations. 

Example  1.  We  compute  the  Euler  equation  with  Lax’s  initial  data.  ut  +  f(u)x  =  0  with  u  = 
(p,  pv,E)T,  f{u)  =  {pv,  pv2  +  p,  v(E  +  p))T,  p  =  (7  —  1  ){E  —  \pv 2),  7  =  1.4.  Initially,  the  density 
p.  momentum  pv  and  total  energy  E  are  0.445,  0.311  and  8.928  in  (0,  0.5);  and  are  0.5,  0  and  1.4275 
in  (0.5, 1).  The  computed  results  by  CO-DG-hrle-3  and  CO-DG-hrlm-3  are  shown  at  T  =  0.16  in 
figure  4,  with  Ax  =  1/200,  A rn  chosen  with  a  CFL  factor  0.4,  A tn  =  0.5Arn.  We  observe  that  the 
resolution  is  quite  good  with  very  small  over/undershoots,  perhaps  among  the  best  within  the  results 
obtained  by  component-wise  limiters.  The  only  concern  is  that  the  contact  discontinuity  is  much 
more  smeared  than  that  of  the  regular  third  order  DG  scheme  with  a  TVB  limiter  (Figure  20  in  [9]). 
We  hope  to  improve  this  performance  by  reducing  the  usage  of  the  reconstruction  limiter  through  a 
troubled-cell  indicator  in  future  work. 

Example  2.  The  Woodward  and  Colella  blast  wave  problem  [43]  for  the  Euler  equation  computed 
by  CO-DG-hrle-3.  Initially,  the  density,  momentum,  total  energy  are  1,0,2500  in  (0,0.1);  1,0,0.025 
in  (0.1, 0.9);  and  1, 0,  250  in  (0.9, 1).  The  density,  velocity  and  pressure  profiles  are  plotted  in  figure  5 
for  T  =  0.01  and  T  =  0.038.  The  solid  line  reference  solutions  are  computed  by  a  3rd  order  central 
scheme  on  overlapping  cells  [31]  on  a  much  refined  mesh  (Ax  =  1/2000).  A rn  is  chosen  with  a  CFL 
factor  0.4,  A tn  =  ^ Arn.  We  observe  stable  results  with  good  resolution  for  this  very  demanding 
problem  in  terms  of  numerical  stability. 

Example  3.  Shu-Osher  problem  [39].  It  is  the  Euler  equation  with  an  initial  data 

(p,v,p)  =  (3.857143,2.629369,10.333333),  for  x  <  -4, 

(p,v,p)  =  (1  +  0.2 sin(5x),  0,  1),  for  x  >  — 4. 

The  density  profiles  are  plotted  at  T  =  1.8,  with  Ax  =  1/40,  see  figure  6.  Arn  is  chosen  with  a  CFL 
factor  0.5,  A tn  =  0.5Arn.  We  observe  very  good  resolution  for  this  example. 
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Figure  4:  Lax’s  Problem.  Ax  =  1/200.  From  left  to  right,  top  to  bottom,  (1)  density  (CO-DG-hrle-3); 
(2)  velocity  (CO-DG-hrle-3);  (3)  pressure  (CO-DG-hrle-3);  (4)  density  (CO-DG-hrlm-3). 


Example  4.  2D  Riemann  problem  [25]  for  the  Euler  equation  computed  by  CO-DG-hrle-3.  The 
2D  Euler  equation  can  be  written  as 

ut  +  f(u)x  +  g(u)y  =  o,  u  =  (p,  pu ,  pv,  E)T,  p  =  (7  -  1)(E  -  \p(u2  +  v2)) 
f(u)  =  (pu,  pu2  +  p,  puv,u(E  +  p))T ,  g(u)  =  (pv,  puv, pv2  +  p,  v(E  +  p))T, 

where  7  =  1.4.  The  computational  domain  is  [0, 1]  x  [0, 1].  The  initial  states  are  constants  within 
each  of  the  4  quadrants.  Counter-clock- wisely  from  the  upper  right  quadrant,  they  are  labeled  as 
(pi,Ui,Vi,pi),  i  =  1,2,  3, 4.  Initially,  p\  =  1.1,  u\  =  0,  v\  =  0,  p\  =  1.1;  p2  =  0.5065,  U2  =  0.8939, 
V2  =  0,  P2  =  0.35;  /?3  =  1.1,  U3  =  0.8939,  V3  =  0.8939,  P3  =  1.1;  and  p 4  =  0.5065,  U4  =  0,  tq  =  0.8939, 
P4  =  0.35.  The  density  and  pressure  profiles  are  plotted  at  T  =  0.25  in  figure  7,  with  100  equally 
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Figure  5:  Woodward  and  Colella  blast  wave  problem  computed  by  CO-DG-hrle-3.  Ax  =  1/400.  Top: 
density;  middle:  velocity;  bottom:  pressure.  Left:  T  =  0.01;  right:  T  =  0.038. 


spaced  contours.  The  numerical  resolution  is  quite  good  for  this  problem.  The  density  profile  along 
y  =  1/3  is  plotted  in  figure  8.  There  is  no  oscillation  near  the  discontinuities. 


Example  5.  Double  Mach  reflection  [43]  computed  by  CO-DG-hrle-3.  A  planar  Mach  10  shock  is 
incident  on  an  oblique  wedge  at  a  7r/3  angle.  The  air  in  front  of  the  shock  has  density  1.4,  pressure 
1  and  velocity  0.  The  boundary  condition  is  described  in  [43].  The  density  and  pressure  profiles  are 
plotted  at  T  =  0.2  in  figure  9,  with  30  equally  spaced  contours.  Ax  =  Ay  =  1/120,  A rn  chosen  with 
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Figure  6:  Shu-Osher  Problem,  Ax  =  1/40.  Left:  CO-DG-drlm-3.  Right:  CO-DG-hrle-3. 


Figure  7:  A  2D  Riemann  problem  [25]  computed  by  CO-DG-hrle-3.  Ax  =  Ay  =  1/200,  Left:  density. 
Right:  pressure. 


a  CFL  factor  0.4,  A tn  =  0.99Arn.  We  can  see  in  the  lower  graph  (the  cross-section  density  profile 
along  y  =  1/3)  that  the  computed  result  is  non-oscillatory. 

6.  Concluding  Remarks  and  a  Plan  for  Future  Work 

In  this  paper  we  have  developed  a  central  discontinuous  Galerkin  method  based  on  staggered  over¬ 
lapping  cells,  with  a  numerical  viscosity  which  stays  bounded  when  the  time  step  size  goes  to  zero. 
Time  discretization  is  via  the  standard  TVD  Runge-Kutta  method.  We  have  also  developed  a  multi¬ 
layer  hierarchical  reconstruction  procedure  and  used  it  as  a  post-processing  limiter  for  our  central 
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Figure  8:  A  2D  Riemann  problem  [25].  Density  profile  along  y  =  1/3. 


Figure  9:  Double  Mach  reflection  computed  by  CO-DG-hrle-3.  Ax  =  Ay  =  1/120.  Top:  density 
contours;  Middle:  pressure  contours;  Bottom:  density  cut  along  the  line  y  =  1/3. 


DG  scheme.  The  limiter  is  able  to  maintain  the  original  order  of  accuracy  and  can  effectively  con¬ 
trol  spurious  oscillations  for  discontinuous  solutions.  In  future  work  we  will  generalize  the  method 
to  convection  diffusion  equations,  improve  the  limiter  by  applying  troubled-cell  indicators,  and  also 
study  further  the  hierarchical  reconstruction  procedure  as  a  limiter  for  the  regular  DG  methods  and 
finite  volume  schemes.  A  stability  analysis  and  error  estimates  for  the  central  DG  scheme  as  well  as 
a  comparison  between  the  regular  DG  and  central  DG  schemes  will  also  be  performed. 
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